Scaling behaviour of trapped bosonic particles in two dimensions at finite temperature 



o 

1—5 

in 



C/2 

I 



I 

o 



en 
> 

o 
m 

o 

cn 

o 



X 



Giacomo Ceccarelli, Christian Torrero 

Dipartimento di Fisica deU'Universita di Pisa and I.N.F.N., 

Sezione di Pisa, Largo Bruno Pontecorvo 2, 1-56127 Pisa, Italy 

(Dated: May 7, 2012) 

In the framework of the trap-size scahng theory, we study the scaling properties of the Bose- 
Hubbard model in two dimensions in the presence of a trapping potential at finite temperature. In 
particular, we provide results for the particle density and the density-density correlator at the Mott 
transitions and within the superfluid phase. For the former quantity, numerical outcomes are also 
extensively compared to Local Density Approximation predictions. 
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I. INTRODUCTION 

Recently, optical lattices have acquired primary experi- 
mental importance since they allow the study of the main 
features of systems of well-localized cold atoms: among 
others, the intriguing interplay of thermal and quantum 
effects in bosonic gases can thus be investigated rather 
precisely, with particular attention usually paid to Mott- 
Hubbard transitions [IH5] . A peculiarity of such experi- 
ments is the confining of particles within a limited region 
of the lattice which is normally achieved by introducing a 
trapping potential. This experimental setup can be mim- 
icked theoretically by the so-called Bose-Hubbard (BH) 
Hamiltonian |17j reading 
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where bl and bi are respectively bosonic creation and 
destruction operators, rii is the local density operator, 
/i is the chemical potential, U is the on-site repulsion 
energy and the sum in the first term is over the nearest- 
neighbor sites of a regular d-dimensional lattice. As for 
the trapping potential V{ri) (being r.^ the distance from 
the center of the trap), a common choice is given by 



Vin) 



VT-^ 



(2) 



being I = J^'Pv~^ the trap size. The exponent p is clearly 
even and will be set to p = 2 in the following. Moreover, 
the energy unit will be fixed by setting J = 1 (thus / = 
1/v) while r and / will be measured in units of the lattice 
spacing a ^ and hence dimensionless. 

In the homogenous case (i.e., with vanishing trap) the 
model undergoes quantum transitions between superfiuid 
and Mott-insulator phases depending on the value of fi. 



The introduction of a trapping potential changes the 
phase diagram [Ttillj: not only a truly diverging cor- 
relation length appears only in the limit I — >■ -l-oo [71 18] 
with /i set to the critical values of the corresponding ho- 
mogeneous system, but also the scaling properties of any 
observable generally acquire an extra dependence on the 
trap size I controlled by the trap exponent 9 given by 



(3) 



A frame to handle this involved scaling is provided by the 
trap-size scaling (TSS) theory [T2[T3]. As a benchmark 
example, at a quantum critical point TSS prescribes the 
free-energy density to scale as 

F{^l, T, I, r) = l-'^^'^+'^T{jif'\ Tf\rl-'^), (4) 

with z the dynamical exponent, v the critical exponent 
controlling how the correlation length diverges, r the dis- 
tance from the middle of the trap, fl = ^ ~ fic, and fj,c 
the critical value of the chemical potential. 

TSS has already been applied to the one-dimensional 
(ID) BH model, both at T = [H] and at finite tem- 
perature |15| : in this paper we extend it to the two- 
dimensional (2D) BH model at finite temperature. In- 
deed, 2D systems are relevant not only from a theoreti- 
cal point of view but have also raised experimental inter- 
est [IIISIIIH]. 

At T = 0, the 2D BH model ^ in the hard-core limit 
(see below) undergoes two phase transitions between su- 
perfluid and Mott insulator at /i = 2 and /i = —2^, while 
at finite temperature it is well known that the model also 
develops a Kosterlitz-Thouless (KT) transition ^IQi 20]. 
In this work we are not going to study the latter but 
rather perform quantum Monte Carlo (QMC) simula- 
tions with chemical potential fixed at /i = —2,0,2 only 
while varying T: This is because our aim is to investigate 
the behavior of the model at the quantum T — critical 



a will be set to 1 from now on. 



^ More precisely, the system is in a Mott phase with {rii) = for 
/i > 2, in a superfluid phase for \fj,\ < 2 and in a Mott phase with 
(rii) = 1 for /^ < — 2. The two transitions share the same critical 
exponents i/ = 1/2 and z =2 |17l . 




FIG. 1: A qualitative sketch of the Kosterlitz-Thouless tran- 
sition (sohd hne) in the fJ, — T plane. 



points. This choice of the parameters should avoid any 
crossings of the KT line, as clear from the qualitative 
diagram in Fig. [ll 

In this framework, the scaling of the particle density 



pin) = {n,) , 
and the density-density correlator 

Gij-i.rj) = (niTij) - {ni){nj) 



(5) 



(6) 



will be studied at fixed trap size I and compared with TSS 
predictions. Besides, we extensively study how numerical 
outcomes for the particle density approach their Local 
Density Approximation (LDA) predictions at the Mott- 
to-superfluid transition with non-zero filling and within 
the superfluid phase. 

Since scaling properties are expected to be universal 
with respect to U, we will work in the hard-core (HC) 
limit U — > +00 where the particle occupation number at 
a generic lattice site can be equal to or 1 only. This 
considerably simplifies the simulation algorithm (which 
is based on the stochastic series expansion [2Tti23] ). 

This paper is organized as follows. In Sec. II we pro- 
vide some details on numerical simulations while in Sec. 
Ill we start our analysis by studying the (rij) =0 Mott 
transition and compare QMC outcomes with the TSS 
theory. In Sec. IV LDA is numerically estimated for the 
2D HC BH model and then applied in Sec. V in consid- 
ering the (rii) = 1 Mott transition. In Sec. VI we analyze 
the superfluid phase by paying particular attention to 
the scaling properties close to those lattice sites where 
the effective chemical potential /ie//, which will be de- 
fined later, equals approximately 2. Finally, we conclude 
in Sec. VII. 



II. QUANTUM MONTE CARLO SIMULATIONS 

Numerical simulations relied on the directed loop al- 
gorithm stemming from the stochastic series expansion 



method [3T]: for a generic system with Hamiltonian H, 
its starting point is given by the standard power series 
expansion of the partition function Z, that is 
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(7) 



a n—0 



being {\a)} a basis set. If H can be decomposed as a 
sum of bond operators Ha^^a- where a^ labels the bond 
and Cj refers to whether the operator is diagonal (q — 1) 
or not (ci = 2) with respect to {\a)} -, Eq. ([t]) can be 
rewritten as 
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(8) 



with Sn standing for a sequence Sn — [ai, Ci], . . . , [a„, c„]. 
We can easily arrange for this setup with the HC BH 
model. The basis {\a)} is chosen to be the set of eigen- 
vectors of the local density operators rii and this auto- 
matically determines which terms in the Hamiltonian are 
diagonal: contributions with blbj have c — 2 while those 
written in terms of the n^'s have c = 1. Moreover, the 
former are already bond-like while the latter have to be 
rewritten: as an example. 



/x ^ n, -> Ai ^ (^ 



(«j> 



Ji Jj 



(9) 



where the sum on the right-hand side runs on nearest- 
neighbor sites and where fi and fj are the number of 
links having respectively site i and site j as one end. 

Even though the Taylor expansion above converges 
|23j , statistically relevant contributions are basically pro- 
vided by configurations where the number of bond opera- 
tors in Eq.(8) is finite and below an opportune value iVtr; 
therefore, truncating the series at order iVtr for practi- 
cal purposes should not entail any significant truncation 
error, as explained in Sec. HA of '22'. In determining 
iVtr, we opted for the standard definition, that is we set 
iVtr = 1-5 Mmax iVbonds/r, whcre M^ax IS the highest ma- 
trix element of the single-bond Hamiltonians and A'bonds 
is the number of interacting site pairs. Besides checking 
that this cutoff was never crossed during the updating 
process, fluctuations in the order of the series expansion 
were monitored to control whether the averaged order 
was consistently less than TVtr (with deviations propor- 
tional to the square root of the mean value) . As proven 
in [22 , these criteria ensure that the truncation error is 
negligible compared to the statistical uncertainty stem- 
ming from Monte Carlo fluctuations. 

Exploiting this truncation, the expression for the par- 
tition function can be further simplified, i.e., 






(-/3)"(iVtr-»)! , 



Ntr 



J]iJ,.,,Ja), (10) 



where A^tr — "- identity operators have been inserted in 
all possible ways in the sequence S^tr- If is understood 




FIG. 2: (Color online) The particle density at /u = 2 with 
T = Tl — 8 and t — 2 for some values of the trap size I. 
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FIG. 3: (Color online) The density-density correlator at fi = 2 
with fixed t = Tl = 8 for different values of the trap size I. 



that now the index Ci can assume a third value (ci = 0) 
corresponding to the identity itself. 

In Eq. ( 10 ) the space of configurations have been gen- 
erahzed to be {|a)} (8 {■S'Art,}. This can be sampled by 
means of two kinds of steps: the first type Mi consists 
of replacing identity operators in the sequence S'^Vt, with 
diagonal ones (and vice versa), while the second kind 
M2 is given by exchanging diagonal operators with non- 
diagonal ones (and vice versa) . An exhaustive description 
of both steps and of how they are performed can be found 
in Ref. [22] ■ Let us just recall that, in implementing kind 
M2, a set of transition probabilities is needed and must be 
determined by solving so-called directed loop equations. 
Depending on the parameters of the Hamiltonian, it is 
possible to select solutions able to reduce the number of 
bounces, '^ that is to cut the amount of moves where the 
proposed change is rejected. Such solutions are preferred 
since they shorten the computer time needed to update 
the configuration. In our simulations, one bounce was al- 
lowed within each group of equations at ^ = and /i = 2 
while two bounces entered into play when /i = — 2. 

One MC step is made out of a single step of type 
Ml followed by a number iVioops of updates of kind M2. 
A'loops is fixed at runtime by imposing that the number 
of visited vertices is of the order of N^^ in a MC step. 

Runs are performed fixing temperature T, chemical po- 
tential /u, trap size I and lattice size L with open bound- 
ary conditions. Finite-size effects are avoided by choosing 
L sufficiently large to obtain L — >■ cxi data within statis- 
tical errors. This condition was fulfilled taking L/l sa 3 
when jjL = 2 and /.t = and L/l « 5 when jjl = —2. 

A standard jackknife was employed to assess errorbars, 
each bin being the mean of 10* MC step measurements. 
Typical statistics of our QMC simulations range from 



2.5 X 10'' MC steps for simulations at ^ = 2 to 7.5 x 10^ 
MC steps for simulations at /i = 0. 



III. THE CRITICAL POINT AT ^ = 2 

We now discuss our QMC results for the density and 
the density-density correlator at the Mott-insulator to 
superfluid transition where {rii) — 0. In analogy with the 
singular part of the free-energy density (|4]), the scaling 
ansatz for the two above-mentioned observables * read: 



p{r) ^l-'''V{fil'\Tl'\rl-'') , 



(11) 
(12) 



where the critical exponents for this transition z/ = 1/2 
and z = 2 have been used. Scaling corrections due to 
irrelevant perturbations in l^^ and possible analytic con- 
tributions have been neglected. After setting d = 2, in- 
troducing the scaling coordinates R = rl~^,R' = r'l^^, 
and considering the system at criticality (so that p, — 0), 
Eqs. (11) and (12 1 can be rewritten as 



rpir)^V{T,R), 
l^OG{r,r')^g{r,R,R'), 



(13) 

(14) 



being r = Tl^^ the scaling variable that controls the 
critical behavior of the system. ^ The meaning of Eqs. 
(131 and (Il4| should be pretty clear: A given observ- 



able rescaled with the proper power of the trap size I 
equals a universal function depending on t, R, R', etc. 
Therefore, data obtained via simulations with values of 
the parameters tuned in such a way to keep the argu- 
ments of the function on the right-hand side of Eqs. ( 13 ) 



^ This issue is treated in great detail in Sec. IID of [22) where the 
XXZ model is studied. Since its matrix elements are in one-to- 
one correspondence with those of the HC BH model, the discus- 
sion can be easily adapted to the present case. 



From the conventions introduced after Eq. ([2|l, it is clear that 
both p{r) and G{r, r') are dimensionless quantities. 
From now on, quantities R, R' and t will always be defined as 
in this section unless differently specified. 



and ( 14 ) constant should collapse on a unique curve once 
that the proper rescaling has been performed. For the 
Mott-insulator to superfluid transition in the low-density 
regime, this condition is fulfilled by performing simula- 
tions with fixed Tl since 9 = 1/2. 

While in the ID HC BH model the particle density 
and the density-density correlator could be treated an- 
alytically both at zero and finite temperature (so that 
numerical outcomes could be compared with their ana- 
lytical values [15]). in the two-dimensional case no exact 
solution is available. In this study TSS is applied to a 
2D-system. 

Figure |2] shows the rescaled particle density. Data in 
it are divided into two groups corresponding to simula- 
tions performed with fixed r = 2 or r = 8. Since t = Tl 
and since the values of I are essentially the same in both 
groups, sets with t = 2 are generally related to lower 
temperatures than those with r = 8. While the latter 
shows scaling corrections at small I, it is evident that the 
former have a more pronounced tendency to collapse on 
a universal curve. This comes with no surprise since uni- 
versality is a feature appearing in proximity of a phase 
transition, which occurs at T = when working with 
chemical potential fixed at /x = 2 as in the present case. 

Figure [3] contains the rescaled density-density correla- 
tor at fixed r = 8 vs. R. In analogy with the particle 
density, also for this observable, numerical outcomes after 
the rescaling prescribed by TSS display a tendency to col- 
lapse on a unique curve when increasing Z, in agreement 
with the ansatz in Eq. (14 1. Once again, corrections can 
be noticed only at smal 



values of the trap size. 



IV. LOCAL DENSITY APPROXIMATION 

In many statistical systems featuring an external po- 
tential V{r) varying with the space position, it is common 
to approximate the ground-state density at point r with 
the value that the density assumes in the homogeneous 
system provided with a constant potential fixed every- 
where at the value V{r) that the potential takes at point 
r itself in the inhomogeneous case. This approximation 
is called Local Density Approximation (LDA). 

LDA has already been verified to be exact in the ID HC 
BH model at zero temperature [T3] and it is reasonable to 
test to which extent it works also in the two-dimensional 
case at finite T. In general, considering a constant poten- 
tial in Eq. (IT]) essentially means to introduce an effective 
chemical potential /icff (j') given by 
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FIG. 4: (Color online) Numerical outcomes for p,{fi) vs. fi 
for different values of the lattice extent L and temperature T. 
Data were collected in the homogenous system with periodic 
boundary conditions. The dotted line represents the polyno- 
mial fit of the data corresponding to the largest extent. 



where p* (/i) is the unknown T = density of the 2D 
homogeneous system provided with an effective chemical 
potential given by Eq. (15). 



In order to obtain an estimate for p* (/i), we performed 
simulations of the 2D system without a trap and with 
periodic boundary conditions^ in the low-temperature 
regime for different values of the effective chemical poten- 
tial. In particular, we employed a set of equally-spaced 
values covering the range from —2 to -1-2. This setup is 
easily obtained by setting to zero the trap parameter v in 
our QMC code and by implementing the specific topol- 
ogy, all other features of the simulation remaining the 
same. 

More specifically, we performed simulations with L = 8 
at T = 1/64, with L = 16 at T = 1/128 and with L = 32 
at T = 1/256, being L the extent of a square lattice. 
Following the same criteria of [24], we checked that the 
data were consistent within errorbars so that we could 
safely assume that the results at T = 1/256 correspond 
effectively to the zero-temperature values. Figure |4] dis- 
plays p* (/i) for the three sets of simulation parameters 
mentioned above; data basically overlap. 

Finally we fitted the T — 1/256 outcomes to a generic 
polynomial function of degree n: 



P* (Ai) = Yl ^' ^' 



(17) 



MGff(?') 



M 



(15) 



Therefore, in analogy with the ID HC BH model, we 
assume that the LDA of the 2D trapped system equals 



where n was chosen by truncating this Taylor expansion 
when the x^ of the fit stabilized. This was the case with 
n = 7 (the reduced x^ being approximately 1.5) though 
truncations at higher order were also considered without 



PLBA{r) 




for fics{r) > 2 , 

for -2 < peff(r) < 2 

for ficff{r) < -2 , 



(16) 



This is expected to reduce finite-size corrections. 



showing meaningful deviations. As expected on theoreti- 
cal grounds, it turned out that the constant term cq read 
1/2 (within 10~^) while even terms were negligible; thus, 
the only non-trivial contributions are given by the odd 
powers for which the following estimates were obtained: 

ci = -0.20779(1) , C3 = -0.01323(1) 

C5 = +0.00441(1) , cj = -0.00093(1) . (18) 



The function in Eq. (17) with n — 1 and coefficients 
as given above is plotted in Fig. |4] and was used for the 
data analysis reported in the following sections. 



n = 1 MOTT TRANSITION 




The invariance under the particle-hole exchange entails 
a similar behavior of the homogeneous HC BH model at 
the transitions with /x = 2 and ^ = — 2. However, the 
trap-size scaling behavior at the (n^) — 1 transition is ex- 
pected to be different than in the vacuum-to-superfluid 
one because the particle-hole symmetry does not hold for 
a trapped system. 

Studies of the trapped ID HC BH model at zero and 
finite temperature [HI [15] have revealed how, at the su- 
perfluid to Mott transition with non-zero filling, the par- 
ticle density approaches its Local Density Approximation 
(LDA) in the large-Z limit. In analogy with the one- 
dimensional case, we thus expect this observable to be 
given by an expression like 






') 



(19) 



since 9 = 1/2 and z = 2 again; irrelevant corrections in 
l~^ have been neglected once more. Therefore, the scaling 
quantity should not be the density itself but rather the 
difference 



Ap{r) = p{r) ~ pi^BAirl ^) 



(20) 



Figure [5] shows how the particle density converges to 
the LDA in the two-dimensional model. Since t = Tl 
as in Sec. HI, once again data sets with t = 2 corre- 
spond to temperatures lower than those of the sets with 
T = 8, given the common values of the trap size. Besides 
improving with increasing I in agreement with Eq. ( |19[ ), 
the convergence to the LDA is better at small T, as clear 
from comparing the upper and middle part of Fig. [5] 
since LDA itself is approached at T — )■ 0. 

The lower part of Fig. [5] illustrates the behavior of 
Ap(r) at T = 8 after the rescaling suggested by Eq. ( 19 ) 



has been performed: a tendency to collapse on a univer- 
sal curve is evident in a region close to the origin while 
some transition-like peaks appear at a distance r fa 21 
from the center, thus drifting with increasing L There- 
fore, in the l^+oo limit, these peaks disappear and only 
the region with the universal curve is left, hence confirm- 
ing TSS predictions in Eq. ( 19 ). The deviations from the 
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FIG. 5: (Color online) The particle density at /i = —2 with 
fixed T = Tl = 2 (top) and t — 8 (middle) for different values 
of the trap size I and the scaling of the subtracted particle 
density at /x = —2 with r = 8 (bottom). The dotted line in 
the first two plots represents the numerical estimates of the 
LDA. 
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FIG. 6: (Color online) Scaling of the density- density correla- 
tor a,t jj, — —2 with fixed r = 8. 



expected scaling at finite I will be treated in Sec. VI. 

As for the density-density correlation function, scaling 
expectations, again given by Eq.( 14 1, are also nicely con- 
firmed, as depicted in Fig. |6] At t = 8, besides scaling 
corrections at small trap size, no strong deviations from 
a universal curve could be detected within errorbars. In 
the ID HC BH a striking result was the universality for 
the correlator between Mott-to-superfluid transitions at 
{ui) — and at (n^) = 1. A closer inspection at Figs. Is] 
and [6] reveals how this feature seems to hold in the 2D 
case as well. 

In the ID HC BH at zero temperature [13], a pecu- 
liar characteristic of scaling quantities like the subtracted 
particle density or the density-density correlator was the 
appearance of modulations depending on the trap size 
I. It was shown that, at fixed /i < 1, there are values 
of the trap size /, whose number increases with I itself, 
for which the energy gap AE between the ground state 
and the first excited state vanishes. This repeated level 
crossing leads to modify the scaling ansatz like Eqs. ( 13 ) 



and ( 14 1 in order to include a further dependence of the 
universal function on a phase (j) related to the difference 
between values of I corresponding to zeros of AE. At 
finite temperature, it is expected that this phenomenon 
plays a lesser and lesser role with increasing T since ther- 
mal fluctuations should prevail on quantum effects. This 
has already been checked in the ID BH model [TS] and 
is confirmed also in two dimensions as Figs. [5] and [6] 
reveal. ^ 
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the rescaled subtracted particle density (bottom) at /i = 
with fixed t = Tl = 2. 



VI. THE SPATIAL REGION WHERE ^eff = 2 

We are now going to study the 2D BH model at fi — 0. 
Since this corresponds to the deep interior of the super- 
fluid phase, no transition is expected to be monitored. 
However, a closer look at Eq. ([I]) suggests to group the 
last two terms of the BH Hamiltonian, thus giving rise 
to the effective chemical potential fj-csir) already intro- 
duced in Eq. (15). 

After this step, it turns out that, in those regions 
where /Xeff(r') « 2, the Hamiltonian is effectively given 



by that of the homogeneous system at criticality; ^ there- 
fore, we might expect some sort of phase transition and 
universal behavior in these regions even if /x is not set to 
a critical value. 

This seems to be indeed confirmed by the behavior of 



^ Recalling the phase diagram, it has to be borne in mind that the 
condition fi < 1 in ID corresponds to ^ < 2 in 2D. 



* The number of dimensions is understood to be d = 2, otherwise 
a critical regime would appear in different sectors of the lattice. 
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FIG. 8: (Color online) The density-density correlator at /j = 
with fixed t = TI = 2. 



the particle density portrayed in Fig. It] If Ap(r) de- 
fined in Eq. (20) is rescaled with I and plotted versus 
r/l as suggested by Eq. (19), ^ no substantial difference 



with the LDA is observed except for some transition-like 
peaks at "critical" distance re = V^l, corresponding to 
the distance from the center obtained from Eq. ( 15 ) after 
setting /Xoff(r) = 2 and ^ = 0. 

We have already encountered such a situation in Sec. 
V. Indeed, by setting ^eff('') = 2 and ^ — —2 again in 
Eq. (15), this results in r = 21, corresponding to the 
distance at which peaks were observed in the {n.i 



= 1 



Mott-insulator phase (see lower part of Fig. [5]). 

In order to study the scaling in such spacial sectors of 
the system, an expansion of //off (r) around Tc is needed. 



^effW = 2+p(2-M)i- 



-i/p 



I 



+ 0[ir~r,f] . (21) 



As pointed out in Ref. [TS] , the length scale ^ of these 
critical modes should behave like ^ f^ l'^ , a being the ex- 
ponent associated to a linear potential. In other words, 
cr = 1/3 |14| and, again neglecting irrelevant contribu- 
tions in l~^ , scaling should read 



Z2/3Ap(r)«P(r,i?), 



(22) 
(23) 



where R and r correspond now to R = {r — rc)/l^'^ and 
T = Tp/^. In deriving the latter exponent the value 
z = 2 has been employed since we expect critical modes 
close to Tc to be controlled by the same dynamical ex- 
ponent as in the regular Mott-to-superfluid transitions. 
The exponents in Eqs. (22) and (23) do not coincide 
with those in Fig. [71 where also T and / are not properly 



® Recall that z = 1 in the superfluid phase while 9 = 1 for smooth 
modes according to an ansatz verified in Ref. |14| . 
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FIG. 9: (Color online) Scaling of the subtracted particle 

density (top) and the density-density correlator (bottom) at 

/i = with fixed Tl^'^ = 1 around distance re where /icff ~ 2. 



tuned, and this explains why data sets do not collapse 
on a unique curve around Vc = ^/2l. Because of the "im- 
proper" values of the trap size and the temperature, also 
the particle-particle correlator G'(r, 0) in Fig. Isldoes not 
show any particular scaling but simply vanishes after a 
few lattice spacings. 

Figure I9] shows the behavior of Ap{r) and G{r,rc) vs. 
(r — rc)/r^^ at TP^^ — 1 after these observables have 
been rescaled according to Eqs. (l22|) and (23). The col- 



lapsing of both quantities on a single curve proves the 
foreseen scaling clearly right. 



VII. CONCLUSIONS 

We have studied the scaling properties of the two- 
dimensional Bose-Hubbard model at finite temperature 
in the presence of a trapping potential at the Mott- 
insulator to supefluid transitions. The interest in this 
system and its properties is not only theoretical given 
that 2D experiments involving cold atoms in optical lat- 



tices are currently being carried out [H [TBI IIH] ■ 

Particular attention has been paid to the particle den- 
sity and the density-density correlator, both computed 
by means of QMC simulations in the hard-core limit 
U — > -foo. The latter choice is motivated by the fact 
that the on-site coupling U should play no role in deter- 
mining the main properties of the system at criticality. A 
comparison between numerical outcomes and TSS ansatz 
has subsequently been performed, revealing that theoret- 
ical expectations are well-motivated. 

An interesting feature arising from our study is how 
LDA compares to the particle density p{r) at finite tem- 
perature in two dimensions. As in the one-dimensional 
case [H], it turns out that also in the 2D HC BH model 
p{r) rapidly converges to the LDA when increasing / even 
at T > 0, regardless the value of the chemical potential 



/i. Even if LDA turned out to be broken in frustrated 
systems |25j and some shortcomings were proven also in 
describing the T = phase diagram of the BH model in 
two dimensions [2Sj, this work shows that LDA is capa- 
ble of describing at least some properties of the 2D BH 
model quite nicely, also outside its theoretical range of 
application. 
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